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Abstract. We present an alternative analysis of CMB time ordered data (TOD) using a wavelet-based representa- 
tion of the data time-frequency plane. We demonstrate that the wavelet transform decorrelates l//-type Gaussian 
stationary noise and permits a simple and functional description of locally stationary processes. In particular, this 
makes possible the generalization of the classical algorithms of map making and CMB power spectrum estimation 
to the case of locally stationary 1/f type noise. As an example, we present a wavelet based algorithm for the 
destriping of CMB-like maps. In addition, we describe a wavelet-based analysis of the Archeops data including 
time-frequency visualization, wavelet destriping and filtering of the TOD. These filtered data was used to produce 
polarized maps of Galactic dust diffuse emission. Pinally, we describe the modeling of the non-stationarity on the 
Archeops noise for the estimation of the CMB power spectrum. 
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1. Introduction 

After their discovery by the COBE satellite 
ISmoot et al. 199^ the Cosmic Microwave Background 
(CMB) temperature anisotropies have become one 
of the most powerful observational tests for cosmol- 
ogy. Indeed, recent measurements of their angular 
temperature power spectrum by ground based ex- 
periments such as DASI IIHalverson et al. 2flfl2|l . CBI 
UPearson et al. 2002|l and VSA ijDickinson et al. 200411 
and by balloon-borne experiments such as BOOMERang 
JNettertiel d et al 200211, MAXIMA HT.ee et nimn^ and 
Archeops IIBenoit et al 200.Sbl ITVistram et «r2005l have 
provided strong constraints on the main cosmological 
parameters such as the Hubble constant, i/o, the Universe 
total density, Oq, matter density, fim, and energy density 
J7a, and the scalar ris, and tensor nt spectral indexes 
associated to the power spectrum of the primordial fluc- 
tuations. More recently, the satellite mission WMAP ^ 
ijBennett et al. 2008|l has measured both the temperature 
CMB angular power spectrum and the cross correlation 
between temperature and polarization l |Kogut et al 2003| l 
leading to high precision determination of the cosmologi- 
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cal parameters. In the near future, the PLANCK satellite 
mission ^ will significantly improve the accuracy on the 
measurement of the angular power spectrum of the CMB 
temperature anisotropies and will also measure the CMB 
polarization anisotropies. This will lead to a rather more 
selective test for discriminating competing cosmological 
models and to constraints on the inflationary paradigm. 



The analysis of the PLANCK satellite data will 
require both faster and more accurate algorithms than 
those currently available, to deal with such large data 
sets containing few 10^ samples, and to measure po- 
larized fluctuation signals which are two or more order 
of magnitude smaller than temperature fluctuations. 
One of the most important problems in the analysis 
of CMB data is the characterization of the properties 
of the noise in the time ordered data (TOD) and 
their exploitation for optimal map-making and for the 
reconstruction of the angular power spectrum of the 
CMB temperature and polarization anisotropies. The 
extremely large size of the data sets makes compu- 
tationally unfeasible the direct inversion of the noise 
covariance matrix and even its storage jBorrill (1999)] |. 
This leads to the common use of Fourier techniques 
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° To be la unched by ESA in 2007 . 
|http: / / sci.esa.int /science-e /www / area / index. cfm?fareaid= 1 7| 



2 



Macias-Perez, J.F. & Bourrachot, A.: Wavelet analysis of time ordered data 



ijDore et al 2nnil I Hi von et al.. Mm lYvon et al pUfB 
which are exclusively valid in the case of stationary 
and Gaussian noise for which the covariance matrix is 
Toeplitz and it is diagonalized by the Fourier transform. 
In the case of PLANCK, observations will take place over 
two years and therefore it seems reasonable to expect 
that the TOD noise properties will change along the 
mission. Further, residuals from badly removed parasitic 
noises will also contribute to the non stationarity and 
non Gaussianity of the TOD noise. 



HBenoit et al 2nn2|l . 

The plan for this paper is the following. Section 2 
introduces the wavelet and wavelet packet transforms. 
Section 3 discusses the properties of the wavelet trans- 
form of stationary Gaussian process. Section 4 deals with 
non stationary noise and in particular locally stationary 
noise. Section 5 discusses a wavelet based destriping al- 
gorithm. Finally, section 6 describes a wavelet analysis of 
the Archeops TOD. 



New promising Fourier methods (|Natoli et al 20011 
based on the concept of locally stationary noise have 
been presented to deal with the non stationarity problem. 
However, they are not well adapted to the case of resid- 
uals from parasitic signals, and they reduce significantly 
the number of data samples to be used in the analysis 
so that the overall instrument sensitivity also decreases. 
In addition, these new methods requires an a priori 
knowledge of the zones of locally stationarity on the data 
but no algorithm has been yet proposed to address this 
issue. In this paper, we describe a new approach to deal 
with non stationarity in CMB data sets. This approach 
is based on the simultaneous analysis of the time and 
the frequency properties of noise using the wavelet 
and wavelet packet transforms. By contrast to Fourier 
based methods, non stationarity is naturally described 
by the wavelet transform and no cutting of the data is 
needed to properly define locally stationary processes. 
In addition, we discuss how wavelet techniques are very 
useful to identify systematics in the TOD both visually 
and automatically. The application of these techniques 
to the Archeops data is presented going from the simple 
visual inspection for classification of the data to specific 
designed algorithms to remove low frequency components 
in the TOD and remaining stripes on the Archeops maps. 

Archeops ^ is a balloon borne bolometer experi- 
ment that aimed at measuring the CMB temperature 
anisotropy over large and small angular scales. It pro- 
vided the first determination of the Ce spectrum from 
the COBE mu ltipoles IjSmoot et al 199^ to the first 
acoustic peak IjBenoit et al 2003al ITristram et al. 200511 
from which it gave a precise determination of various 
cosmological parameters, such as the total density of the 
Universe and the baryon fraction ijBenoit et al 2nn3hjl . 
Archeops was also designed as a test bed for Planck-HFI 
and therefore shared the same technological design: a 
Gregorian off-axis telescope with a 1.5 m primary mirror, 
bolometers operating at 143, 217, 353 and 545 GHz 
which are cooled down at 100 mK by a '^He/^He dilution 
designed to work at zero gravity and similar scanning 
strategy. Because of this Archeops data are expected to 
have common features with Planck data with respect to 
random noise and systematics. A detailed description 
of the instrument and its performances can be found in 



2. The wavelet and wave let- packet transforms 

In this section we give a brief overview of the wavelet and 
wavelet packet transforms which will be extensively used 
in the following sections. 



2.1. The continuous wavelet transform 

The term wavelet is used to refer to a set of orthonor- 
mal basis functions generated by dilation and translation, 
'4'a,b{t) = a~ 2 ip (^^) of a compactly supported function, 
ip{t) where a and b are called the scaling and transla- 
tion factors. ip(t) is called the wavelet function or mother 
wavelet and it is assumed to satisfy the admissibility con- 
dition 

_ _ . 2 



Co! 



dw < 



(1) 



where ip{w) is the Fourier transform of tp{t) 

The continuous wavelet transform (CWT) of the function 

g{t) is defined by 



-t-oo 



Wa,b = {f,i>a,b) 



9{t) Ipa.bit) dt. 



(2) 



where Wa^b are called the wavelet coefficients associated to 
the function g{t). We can also define the inverse wavelet 
transform so that 



+ 00 +00 
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(3) 
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The CWT constitutes a powerful tool to describe simul- 
taneously the time and frequency properties of a given 
function g{t). In this sense, it supersedes the standard 
Fourier transform within the limits of the Heisenberg's 
uncertainty principle A/ At > ^ where / = 2ttw. 



2.2. The Discrete Wavelet Transform 



|http : //www. archeops . org| 



The discrete 
hereafter, of 

{X{tm), m = 
the CWT, Wab 



wavelet transform, denoted DWT 
a discretely sampled time series, 
0, . . . , — 1} can be obtained from 
, setting a = and b = 2^~^ kj for 
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j ~ 1,J and kj 0, . . . , 2^^^ where TV = 2-^™"^. In other 
words, the DWT is obtained by restricting the CWT to a 
'dyadic partition' in time. For each scale Tj = 2^~^, with 
j = 1, . . . , J, the DWT is composed of Nj — |j wavelet 
coefficients W^-. 

Therefore, the DWT is an orthonormal trans- 
form W that takes a time series - TOD - X = 
{Xo,Xi Xn-i}'^ and yields a vector of N DWT 

coefficients 

W = WX= {Wf,W^,...,Wj,Vj}. (4) 

As above, the sub-vector Wj contains Nj = |j wavelet co- 
efficients associated with scale r, = 2^~^, whereas Vj con- 
tains ^ scaling coefficients associated with scale Aj = 2'^ . 
These scaling coefficients are needed to compensate for the 
truncation of the DWT at level J. As it will be discussed 
later on the section, they account for the smooth compo- 
nent of the time series which is left over by the wavelet 
coefficients. Indeed, for J = Jmax = lg2 ^ then Vj^^^ is 
composed by a single scaling coefficient which represents 
the mean value of the time series. The orthonormality con- 
dition yV^W = In impHes that the inverse DWT is given 

by 

X = W^W = {Xo, Xi , . . . , XN-if (5) 

We can formally describe the DWT in terms of 
wavelet and scaling filters as follows. Assume hi ~ 
{hifi, . . . , hix-i, 0, . . . , 0}^ to be a vector of length N 
whose first L < N elements are the unit wavelet 
filter coefficients for a given compactly supported 
wavelet l |Daubechies,I. 1992| | and with discrete Fourier 
transform (DFT) {Hi^k, k = 0, . . . , N - 1}. Then gi = 
{5i_o, . . . 0, . . . , 0}"^ is a vector of length N con- 

taining the zero padded scaHng filter coefficients for unit 
level, defined via 51^; = (—1)'^^ for ^ = 0, . . . , L— 

1 with DFT {Gi^k',k = 0, . . . , TV - 1}. To obtain the level 
1 wavelet and scaling coefficients, X is filtered using hi 
and gi, and then resampled by 2 

min(L,7V)-l 

Wl^n = J2 ^l.l ^2{n+l)~l-l mod n 

min(L,Af)-l 

Vl,n = E 91.1 X2 (n+1) — 1 — / mod n 

/=0 

where n ~ 0, . . . , N/2 — 1 and hi and gi can be 
considered as high pass and low pass filters respectively. 
Therefore, the wavelet coefficients Wi represents the 
coarse component of X and the scaling coefficients Vi 
the smooth part of X . 

To obtain the level 2 wavelet W2 and scaling V2 coef- 
ficients Vi is high pass and low pass filtered as before 
using hi and gx and then subsampled by 2. Finally, at 
level j, Wj and Vj are given by subsampling by 2 the 
high pass and low pass filtering of Vj_i using hi and 
gi. Here, the time series is considered to be circular 
and as a consequence the first L — 2 coefficients at each 
level j are affected by boundary conditions and therefore 



their statistical properties differ from those unaffected by 
circularity. This algorithm is known as the DWT pyramid 
algorithm l |Mallat, 1989| |. 

It is convenient to define equivalent hj and gj filters 
with Lj ~ [2^ — — 1) + 1 non zero elements such that 

min(Lj ,JV)-1 

= hj.l X2i(n+l)-l-l mod n: n = 0, . . . , iVj 

(=0 
min(Lj,N)-l 

^j.n — S 9j.J X2i{n+l)-l-l mod m U = 0, . . . , Nj 

1=0 

(6) 

and which are respectively the inverse DFT of 

Hj^k = Hi2i--^k mod N W ^121^- mod JV > k = Q,. . . , N —1 
1=0 

(7) 

and 

i-i 

Gj,k = Yl_'~^l,2'k mod N, k = 0,...,N-l (8) 

/=o 

This equivalent filter representation permits a clear 
understanding of the overall filtering applied to the time 
series to obtain the wavelet and scaling coefficients at 
level j. Within this framework it is obvious that the 
coarser components of X are represented by the low j 
wavelet coefficients and the smoother ones by the large j 
wavelet coefficients and the scaling coefficients. Indeed, 
the DWT leads to a dyadic sampling in frequency so that 
the wavelet coefficients are obtained by band pass filtering 
of X in the frequency interval lj = ^^rrr — I /I — ^ for 
unity sampling frequency. 

Finally, the DWT can be also represented as the de- 
composition of the time series into an orthonormal base 

J N/2^-1 N/2-'-l 
j=l k=0 k=0 

(9) 

and therefore, the total energy is conserved 
J 

l|X||' = ^||W,||V||V^||2 (10) 

where ^pj^k and (j>j^k are respectively the wavelet and scal- 
ing functions. 

2.3. Wavelet bases 

In this paper we concentrate on compactly sup- 
ported orthogonal wavelet bases such as the Haar, 
Daubechies Coifffets, Least Asymmetric (also called 
symmlets) and Best LocaHzed wavelets which are 
extensively used in the statistical community (see 



4 



Macias-Perez, J.F. & Bourrachot, A.: Wavelet analysis of time ordered data 



for example |Daubechies,I. 1992| IPercival fc Walden 20001 
|Vidakovic,B. 1999[ . The Haar wavelet and scaling filters 
are given by 



J^Haar ^ 1 1/^2, - l/\/2} 

^ {l/V2,l/%/2} 



9 



(11) 



From this definition, it is obvious that the h^"-"-'^ filter is 
a high-pass filter and the wavelet coefficients are given by 
the difference between adjacent data samples. In addition, 
the g^°-°-^ filter is clearly a low-pass filter and the scaHng 
coefficients are given by the mean value of adjacent data 
samples. 

The so called Daubechies, CoifHets, Least Asymmetric 
and Best Localized wavelet bases correspond each of them 
to a different wavelet family. Each wavelet family is de- 
fined to fulfill a particular requirement like for example 
the symmetry of the wavelet function in the case of the 
Least Asymmetric family. For each family the width, L, of 
the scaling and wavelet filters is a variable parameter and 
as discussed by |Lai, M-J. 1995l it is directly related to the 
shape of the low-pass and high-pass equivalent filters. In 
particular, the larger L the closer the low and high pass 
filters will be to a perfect low pass and high pass filter 
respectively. This will be discussed in more details in the 
following sections. 

2.4. The Maximal Overlap DWT 

The Maximal Overlap DWT (MODWT) also known 
as the non-decimated DWT or the stationary DWT, is 
a particularly interesting generalization of the DWT. 
Actually, it allows us to keep the original sampling in 
time of the time series at each level j of the wavelet 
transform. In addition, the MODWT is invariant with 
respect to time shifting while the DWT is not. These two 
properties will reveal important when studying the time 
evolution properties of the timeline. 

From a mathematical point of view the level Jo 
MODWT of a time series Xt of data samples is a 
transform consisting of Jo + 1 vectors Wi, . . . , Wj„ and 
Vj„, all of which have dimension N. The vector Wj 
contains the MODWT wavelet coefficients^ associated 
with changes on scale Tj = 2^~^ while the Vj^, contains 
the MODWT scaling coefficients associated with averages 
on scales Xjg = 2'-'°. 

As for the DWT Wj and Vj can be calculated by 
filtering Xt, namely. 



= Y^^=0 ^ hj^iXt-l mod N 
= X]i=0 9j,l^t-l mod N 



(12) 



t = 0, ...,iV — 1, where and are the jth 

level MODWT wavelet and scaHng filters. These filters 
are defined in terms of the DWT wavelet and scaling 



filters as hjj = hjj/2'i and gj,i = gj.i/2i and have 
width Lj = (2-' — — 1) + 1. In practice the wavelet 
and scaling coefficients are not calculated from equation 
[T2I but via a pyramidal algorithm similar to the DWT 
(IPercival Walden 20001 |Vidakovic,B. 1999| |. Note that 
the first min{Lj — 2, TV — 1} elements in both Wj and 
Vj are affected by circularity and therefore are boundary 
coefficients. 

UnHke the DWT, the MODWT is not an orthonor- 
mal transform because at each level j the components of 
both Wj and Vj are not independent. Nevertheless, the 
MODWT is capable of producing a multi-resolution anal- 
ysis of the data such that 



+ 



v„ 



(13) 



Finally, the MODWTj^n also_^be expressed as^a matrix 
operations such that Wj = WjX and Vj 
time series X can be recovered from 



VjX. The 



X^^WjW.+VjVj 

2.5. The Discrete Wavelet Packet Transform 

The wavelet transform as defined above is just a repre- 
sentation of temporal signals at different scales j which 
correspond to the frequency intervals {^jTt,^} • This 
representation although extremely useful and powerful, is 
unusual from the Fourier point of view. At this respect, 
it is interesting to define the discrete wavelet packet 
transform, DWPT, by generahzing the DWT, so that we 
can obtain an equipartition of the frequency domain and 
of the time domain similar to that given by a Windowed 
Fourier Transform (WFT). In simple words, the DWPT 
is obtained at each level j by high pass filtering with 
{hij} and low pass with {51,;} the smooth component 
Vj_i in the DWT, but also filtering high pass and low 
pass the coarse component Wj_i. 

The level j DWPT of a iV = 2^^ dimensional vector 
X is an orthonormal transform yielding an N dimensional 
vector of coefficients that can be partitioned as 



0,...,2-'"-l} 



where each W, „ has dimension N.i 



2-^-1 and 



IS 



nominally associated with the frequency interval Ij,„ = 
{2^^Tr}- Note that we define Wo,o = X following pre- 
vious definitions. Together these 2^ vectors divide the fre- 
quency interval [O, ^] into 2^ intervals of equal width ^ttt- 
This leads to an equipartition both in time and frequency 
of the time-frequency plane similar to that produced by 
the WFT without explicit time partitioning. 
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Defining Wj^n,t the th element of „ we can write 

L-l 

^3,n,t — 5Z ^J-l,LfJ.2t+l-i mod JVj.i I t ~ . . . ,Nj 

1=0 

(14) 

where 

^ J gi,;,if n mod 4 = or 3; , , 

- \ /ii,,,if n mod 4=1 or 2; ^^^^ 

and J is f if n is even and otherwise it is 

Since the DWPT is an orthonormal transform, we can 
use it to partition the energy in X via 

2^-1 
n=0 

where || Wj „||^ can be interpreted as the contribution 
to the energy due to frequencies in the band X,,n. As the 
above relationship is valid for any given level j in the inter- 
val 1 < j < J, we can consider the DWPT for 1 < j < J 
is an overcomplete representation of X. Further the ba- 
sis functions associated to it constitute a sort dictionary 
which is commonly used for representing functions via 
Matching Pursuit techniques ( |Mallat fc Zhang, Z., 1993| | . 



which using the DFT transforms into 

^ cov{Wj.t W,,^t'} = 

(19) 

where Sx{f) is the spectral density function of Xt. 
Under the assumption of nominal band pass filters, Hj{f) 
and Hj'{f) have band-passes given by Ij and Tj/ which 
for j ^ j' do not overlap and thus the correlation of 
wavelet coefficients across scales should be negligible for 
a Gaussian stationary process. 

Further, when j = j' and t' ~ t + t the correlation 
between wavelet coefficients reads 

Cw{j,r)^E{W,,tW,,t+r}^ 

jK e^^'^'^f^n.U) Sxif) df ^^^^ 

2 

Once again under the approximation of nominal band- 
pass filters we can argue that the above should be 
approximately zero for r ^ when Sxif) approx- 
imately constant in the interval X,-. In other words, 
at each level j the wavelet coefficients Wj behaves as 
uncorrelated white noise the same way Fourier coefficients 
do. 



3. Wavelet analysis of Gaussian stationary noise 

The analysis of stationary Gaussian processes is tradition- 
ally performed in the Fourier domain because of the decor- 
relation properties of the Fourier transform. Indeed, the 
covariance matrix, Cx, of stationary Gaussian distributed 
noise, Xt, is diagonal in the Fourier space and completely 
described by the noise power spectrum, Px, 



Cx{fJ') = diag{Px{f)} 



(16) 



This property makes trivial the calculation of the inverse 
of the covariance matrix and from the point of view of 
CMB analysis, it eases considerably tasks like map mak- 
ing and angular power spectrum estimation. However, 
this feature can not be extrapolated to non stationary 
Gaussian noise. In this section, we show how the DWT 
can also decorrelate Gaussian stationary noise. Therefore, 
based on the localization properties of the DWT we can 
reasonably expect it will also decorrelate non stationary 
Gaussian noise as discussed in the following sections. 

3.1. Correlation between wavelet coefficients. 

Using equation El and assuming a stationary process Xt 
with autocovariance 



SX,r^ E{Xt,Xt,r} 



(17) 



which depends only on r we can write the correlation be- 
tween non boundary wavelet coefficients as follows 

^Lj-l^i,'-l , , (18) 
2^1=0 2^1' '^j-l'^j'-l' ^X.,2i(t+l)-l-2i' (t' + l)+l' 



As discussed above, the transfer function of the wavelet 
filter Tij = Hj H* at each level j can be well approximated 
by a nominal band-pass filter taking the form 



n^if) = 



2J 




< 



2J + - 

otherwise 



(21) 



Within a particular wavelet family, this approximation 
will be more accurate for large filter widths L, however 
the number of boundary coefficients = L — 2 will be 
more important. In practice, it is always possible to find 
a relative low value of L for which the wavelet filter acts 
nearly as a perfect high pass on the data such that the 
decorrelation between wavelet coefficients takes place. 

In summary, for a stationary Gaussian process Xt with 
spectral density function Sxif) ^^^d under the nominal 
band-pass filter approximation the wavelet coefficients at 
scale j are uncorrelated with variance 



a- 



r{W,}=Ji,n,{f) Sx{f)df 

2 

■J^Sxif)df 

2J + 1 



(22) 



We have ignored the scaling coefficients. However, we can 
also calculate their variance from the r.m.s. of the timeline 
so that 



Cj+i = var {Vj} = Nx \ var {Xt} - ^ 



var {Wj } 



23 



(23) 
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Fig. 1. From left to right and from top to bottom, input power spectrum of simulated l//-type noise, histogram of 
the simulated noise wavelet coefRcients and their expected values at scale j = I, correlation of the wavelet coefficients 
and their expected values for scale j = 1 and j = 10 and cross correlation of the wavelet coefficients at j = 1 and 
j = 10 with scales j ~ 2 and j = 11 respectively. 



3.2. Wavelet description of 1/ f-type noise 

For most CMB data sets the noise associated to each de- 
tector timeline corresponds generally to a stationary long 
memory random process also known as l//-type noise. 
This kind of processes present long term correlations in 
time which show up in the Fourier domain as a strong 
increase of power with decreasing frequency. The specific 
characteristics of the noise depends very much on the 
detector itself as well as on the operational conditions 
Hke for example the temperature of the focal plane bath, 
the current and intensity applied to the detector, the 
electronics, etc. In general, the noise properties of given 
detector are badly known and directly derived from the 
data themselves. For large data sets as would be the 
case for the PLANCK satellite surveyor, this task can 
be extremely expensive in computation and requires fast 
techniques and algorithms. In this section we show how 
the wavelet transform can be used to represent and study 
l//-type noise based on the previous results. The DWT 
has two main advantages with respect to the DFT. First, 
from the computationally point of view it is faster and 
presents no fimits in the number of samples to be used as 
the DFT do; and secondly, the wavelet representation of 
stationary process is more compact in terms of number of 
coefficients. 

The simplest 1 / /-type noise is the white noise which is 
completely uncorrelated in the time domain and presents 
a flat spectral density function of the form 



(24) 



Under the hypothesis of nominal band-pass wavelet fll- 
ters and using equation [23 the variance of the white noise 
wavelet coefficients at scale j is given by 



C, 



WN 



(25) 



Using equation [201 we can show that the white noise 
wavelet coefficients are fully uncorrelated at each scale j 



ihr) 



(jWN 

0^ 



otherwise 



(26) 



As the white noise wavelet coefficients are also uncorre- 
lated between different scales j and j' we can state that 
the wavelet transform of a white noise is a also a white 
noise of the same spectral density function. 

In general the spectral density function of l//-type 
noise is represented as 



fk 



(27) 



where fknee is known as the knee frequency. For / > fknee 
the spectral density function is dominated by the white 
noise term which is parametrized by a. By contrast, for 
/ < fknee the 1// dominates. From equation [22 the vari- 
ance of the wavelet coefficients at level j for 1 / /-type noise 
is given by 



[l + 2^+i/fe„eeln(2)] 



1 



if a^l 
if a > 1 
(28) 

Further using equation[2Slthe correlation between wavelet 
coefficients is given by 



cos (2J+i^/t) /-" df 



(29) 



which is largely dominated by the white noise term 
so that the wavelet coefRcients are quasi-uncorrelated. 
This equation is only valid for nominal band-pass filters. 
However, from a practical point of view increasing the 
number of non-zero coefRcients in the wavelet filter will 
flatten the 1/f term in the integral leading to a desired 
level of decorrelation between wavelet coefficients. 

To test the above statements we have simulated a TOD 
of l//-type noise with 2^^ samples. Figure shows from 
left to right and from top to bottom, the input power 
spectrum of simulated l//-type noise, the histogram of 
the simulated noise wavelet coefficients at scale j = 1 , the 
correlation of the wavelet coefficients and their expected 
values for scale j = 1 and j = 10 and the cross correla- 
tion of the wavelet coefficients and their expected values 
at j = 1 and j = 10 with scales j = 2 and j = 11 re- 
spectively. We observe that the wavelet coefficients show 
a Gaussian structure. Further, they are quasi-decorrelated 
as expected within the same scale and totally uncorrelated 
between scales. 



3.3. Simulation of Gaussian stationary noise using 
wavelets. 

From the above discussion we have conclude that un- 
der the nominal band-pass filter approximation the DWT 
decorrelates Gaussian stationary processes. Thus it is pos- 
sible to simulate Gaussian stationary processes via its 
DWT UPercival et al. 20001 [Vidakovic,B. T999| l, the same 
way it is performed using the DFT. Indeed we can easily 
simulate stationary and Gaussian time series Xt of length 
N = 2'' and spectral density function Sx{f) performing 
the following steps 
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Fig. 2. Left panel: simulated l//-type noise using a wavelet algorithm. Right panel: power spectrum of the wavelet 
simulated 1/ /-type noise compared to the input power spectrum in red. 



1. Using a random-number generator, generate a vector 
Zm containing M ~ AN Gaussian white noise deviates 
with zero mean and unit variance 

2. Integrate numerically equation [22l to calculate the 
approximate wavelet coefficient variances C,-, j = 
1, ... ,J+2 and use equation |22l to calculate the 
scaling coefficient variances. 

3. Multiply the first M/2 elements of Zm by ^/CT; the 
next M/4 values by y/C2 and so forth, until to the final 
four elements, which need to be multiply by y^Cj+i, 
\/Cj+i, \/ Cj+2 , \/Cj+3 respectively. 

4. Compute the inverse DWT of the modified Zm vector 
and then randomly choose M samples from it. 

The above recipe can be easily understood in terms 
of the wavelet covariance matrix of a Gaussian station- 
ary time series with N = 2'' samples which due to the 
decorrelation property is given by 

Sw = rfjflff { C'l, . . . , Ci, C2, . . . , C2, ... , 



Cj,. . 




^ of these ^ of these 

,Cj^^i^^j^,C,j,Cj+i } 

2 of these 



(30) 



The left panel of figure |2l shows a simulation of 1 //- 
type noise obtained from the algorithm described above. 
In the right panel, we represent the power spectrum of 
the simulated data which is in agreement with the input 
power spectrum plotted in red. 

4. Locally stationary Gaussian noise. 

Locally stationary processes appear in many physical sys- 
tems in which the mechanisms that produce random fiuc- 
tuations change slowly in time. Over given time inter- 
vals, such processes can be approximated by a stationary 
one. This is the case, for example, of the variations ob- 
served in the noise total power of some CMB detectors 
due to fiuctuations in the focal plane temperature as dis- 
cussed in subsection From the Fourier point of view 
we can imagine the variation of the data power spectrum 
with time. Actually, we could also define locally-stationary 
Gaussian processes as those for which the time-frequency 
plane is divided into time intervals corresponding each of 
them to a stationary process and which are uncorrelated 
between them. Notice that the above definitions of locally- 
stationary processes are very general and in the limit of 
infinitesimal time intervals they lead to non-stationarity 
in a more general manner. 

4.1. Fourier approximation to locally-stationary process 

For a zero-mean random process Xt we can define its time 
evolving autocovariance as 

s^{t,t')^E{X{t) , X{t')} (31) 



which can also be expressed in terms of the distance be- 
tween t and t' and the mid-point between them 



(32) 



Note that under this definition the covariance of a station- 
ary process reads 

cC-^,t-t') = Cit-t')^C{T) (33) 

For a locally stationary process, we expect that in the 
neighborhood of any x £ TZ, there exists an interval of 
size l{x) where the process can be approximated by a sta- 



tionary one. Therefore for any G 

,t + t' 



_ l(x)_ 

X 2 7^1 2 



l(x) 



C{- 



-,t-t') ^Cix,t^t') 



(34) 



Further, if t and t' are far apart and do not belong to 
the same interval of stationarity then we expect that 
X{t) and X{t') are uncorrelated. This means for any 



u e 



if > %i then 



Ciu,v)^E[x{u+^),X{u-^)}^0 (35) 

In the time frequency-plane we can define a time- 
varying spectrum A(t, w) which would be given by the 
Fourier transform of the time-varying covariance C{t + 
|, t — |) with respect to v. 

For stationary processes A{t, w) does not depend on time 
and corresponds to the standard power spectrum. For 
locally-stationary processes ijMallat et al. IQflSjl the spec- 
trum varies with time but it is constant in time within the 
stationarity intervals defined above. This means that the 
same way the Fourier transform decorrelates Gaussian sta- 
tionary processes, the windowed Fourier transform decor- 
relates locally stationary processes. In other words, the 
covariance matrix of locally-stationary Gaussian noise is 
diagonal in an orthogonal base given by performing a 
windowed Fourier transforms in each stationary interval 
l|Mallat e.t al. 1998)l . This base divides the time-frequency 
plane in rectangles of the form 



2 



X + 



l{x) 



where fk = k 



h 



1,2,, 



1 



2 l{x) 



11 



Jk + 



1 



2 l{x) 



2 l{x) 



and 



the sampling frequency. Note that the DWPT produces 
a similar decomposition of the time-frequency plane 
and can also be used to obtain a diagonal base for 
the covariance matrix of stationary Gaussian process 
Xt. For this purpose, best basis algorithms combin- 
ing the DWPT at all available scales j have been 
developed l |Coifman, R.R. fc Wickerhauser, M.V. 19921 
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Fig. 3. From left to right, locally stationary Gaussian noise obtained from two white noises of variances 1 and 25, and 
its wavelet transform with wavelet coefficients ordered from left to right for decreasing scale. 



IPercival & Walden 20001 |Vidakovic,B. 1999| |. We will not 
account for them in this paper. 

In the previous section we have shown that the DWT 
decorrelates Gaussian stationary processes under the ap- 
proximation of nominal band-pass wavelet and scaling fil- 
ters. Therefore, the DWT will also decorrelate a locally- 
stationary Gaussian process which can be considered as 
series of attached stationary Gaussian processes uncorre- 
lated between them. The same way we have defined above 
a time varying power spectrum we can also consider the 
evolution in time of the variance of the wavelet coeffi- 
cients at each level j. As wavelet coefficients are local- 
ized in time, this evolution can be naturally computed 
from the DWT of the data. Thus, for Gaussian locally- 
stationary processes we expect the variance of the wavelet 
coefficients to change from one stationary interval to an- 
other and to be constant within a given interval. This 
statement can be fully understood considering eauationl22l 
and accounting for the time evolution of the data power 
spectrum. To clarify this issue we can write, from equa- 
tion EOl the covariance matrix of the wavelet coefficients 
of a Gaussian locally-stationary process composed of two 
stationary time intervals 



Ew = diag { Cj, C| 



C2 /^i2 
7 1 ^-9 




of these 



C2 y— (2 f^l y— fl 

J- of these ^ of these 



of these 
1 of these 1 of these 




(36) 

where the indexes 1 and 2 refer to each of the station- 
ary time intervals respectively. In summary, the DWT 
decorrelates Gaussian locally-stationary processes, and 
naturally permits the identification of the distinct inter- 
vals of stationarity in the data using the time evolution 
of the wavelet coefficient variance. 

For illustration, we show in the left plot of figure El an 
example of a locally stationary noise given by the super- 
position in time of two white noises of variances 1 and 
25. Its wavelet transform is traced in the right plot with 
the wavelet coefficients ordered from left to right for de- 
creasing scale. We can clearly observe for each scale two 
distinct parts corresponding to the two white noises. 

4.2. Time modulated-stationary processes 

Often, as discussed in the introduction to this section, the 
noise total power of CMB instrument detectors changes 
with time due to slow fiuctuations on the temperature of 



the focal plane bath or drifts in the background power. 
In general, this leads to locally-stationary Gaussian 
noise for which the power spectrum varies along the 
observation period but in the same way for all frequency 
components. Indeed, the total noise power increases 
and/or decreases with time but the shape of the power 
spectrum remains the same all over the observation 
period. These processes can be well approximated 
by a stationary Gaussian process modulated on time 
dFryzlewicz, P., Van Bellegem, S. fc von Sachs, R. 2002| | 
which is defined by 



Xt=a{t) X Yt 



(37) 



where It is a Gaussian stationary process and <7{t) is a 
function of time. From a more physical point of view, 
cr(t) represents the change on the noise total power. For 
convenience we will consider in the following that the 
variance function cr^(t) is defined on (0, 1) by assuming Yt 
has as power spectrum the time-constant power spectrum 
of Xt. Finally, notice that time modulated processes are 
just a particular case of locally-stationary processes. 

From a practical point of view, the above model is quite 
handy because only a(t) has to be estimated from the 
data. Simple although biased estimates of a{t) can be ob- 
tained by comparing the non-time dependent power spec- 
trum of the time series Px{w) to power spectra, Px{w,tk) 
calculated at different but connected time intervals along 
the observation period 



(38) 



This estimate leads to the approximation of a{t) by 
a step-function so that the time intervals have to be 
particularly chosen to represent as best as possible 
the shape of crit). By contrast, using the MODWT it 
is possible to easily obtain unbiased estimates of a{t) 
( |Fryzlewicz, P., Van Bellegem, S. &: von Sachs, R. 2002| | 
from 

J 



W. 



(39) 



The extra factor 2~^ can be easily understood from the 
normalization factor introduced in the definition of the 
MODWT wavelet and scaling filters when deduced from 
the DWT ones. Few extra comments on equation EHI are 
needed. First, a{t) is not defined on (0, 1) and has to be 
renormalized to be compared to the Fourier estimate. 
Second, the MODWT estimate is not consistent and has 
to be smeared out to reduce its variance. Third, we could 
also obtain an unbiased estimator from the DWT in a 
similar manner but its construction is more complex due 
to the time undersampling performed by the DWT as 
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increasing levels. 

In summary, the time-modulated stationary model for 
random processes can be of great interest for dealing 
with Gaussian locally-stationary CMB time series which 
present time evolution in the noise total power but not in 
their power spectrum shape. Indeed, once estimates for the 
a{t) function are available the simulation of such processes 
is reduced to the simulation of a stationary Gaussian pro- 
cess from the power spectrum of Xt which is then multi- 
plied by that function. Further, the covariance matrix of 
such processes can be well approximated in the wavelet 
space as shown in the following subsection. 

4.3. Time modulated-stationary wavelet processes 

Based on the definition of time modulated-stationary pro- 
cesses and on the properties of the wavelet coefficients of 
stationary and locally stationary processes we can also de- 
fine time modulated-stationary wavelet processes as those 
for which the variance of the wavelets coefficients at level 
j is of the form 

var{Wj^t} = a^{t) xCj (40) 

where Cj is a constant and aj{t) is a function of 
time. Notice that this definition accounts also for time 
modulated-stationary processes for which the different 
<Tj{t) functions are a consequence of the undersampling 
performed in the DWT. 

As in the above definition we have supposed the 
wavelet coefficients are uncorrelated their covariance ma- 
trix is therefore diagonal of the form 

Sw = diag { o-i(t}) Ci,ai{tl) Ci, . . . , 
V ' 

of these 

^ . ' ^ , ' (41) 

^ of these -~ of these 

2 of these 

where represents the fcth time sample of the level j 
DWT. 

4.4. Statistical analysis of CMB data sets with locally 
stationary noise. 

For most of the time domain processing of CMB data sets 
Hke for example subtraction of systematics and optimal 
map making we have to deal with the inversion of the 
noise covariance matrix. In general, we consider the noise 
Gaussian and stationary or piece-wise stationary, so that 
the noise correlation matrix is diagonal in Fourier space. 
Therefore, this can be trivially inverted for each piece of 
data. Above, we have shown that the covariance matrix of 
correlated Gaussian and non stationary time modulated 



noise (in particular locally stationary) is diagonal in 
the wavelet space and can be simply described by a set 
of coefficients Cj, where j is the wavelet scale index, 
modulated by a time dependent function aj{t). Therefore, 
the inversion of this matrix is also trivial in wavelet space. 

If we consider Maximum Likelihood algorithms like for 
example optimal map making, the resolution of the sys- 
tem involves terms of the form N~^d where N is the noise 
correlation matrix and d a data vector. In Fourier space we 
perform this operation by first taking the Fourier trans- 
form of d, then dividing it by the diagonal terms of the 
Fourier representation of N and finally transforming back 
to real space. Equivalently, in the wavelet case we can 
write for the stationary case 

N-'d = W^ [{Wdj^ JC,) (42) 

where W represents the direct wavelet transform, the 
inverse wavelet transform and Cj are the variances of the 
wavelet coefficients of the noise wavelet transform. In the 
same way for the non stationary time modulated noise we 
can write 

N-^d = [{W (d/a(t))}^,, /(Q)) (43) 
where cr{t) is the time modulation function. 

In general the noise covariance matrix is not known 
and has to be estimated directly from the data. For the 
stationary case, this is simple as we only need to compute 
the variance of the wavelet coefficients for each of the 
different scales. For non stationary time modulated noise 
we can have a very good approximation using Eq. 1301 by 
first computing the variance of the wavelet coefficients 
assuming stationary noise and then computing the time 
modulation function from the global evolution of the 
wavelet coefficients. 

The top panel of figure, 21 shows in black the global 
time variation of the noise of one of the Archeops bolome- 
ters which decreases with time (a complete description 
of the Archeos data is presented in Section ^ . In red we 
overplot the time modulation function, a{t) as determined 
from the data themselves following the above technique. 
Finally, the bottom panel of the figure shows the global 
time variation of a simulation of the Archeops noise. This 
have been obtained by multiplying by a{t) a realization 
of a stationary Gaussian noise produced as described in 
section I301 In red we overplot the a(t) estimated for the 
simulation which behaves as the one estimated from the 
data. Notice that this kind of non stationary simulations 
of the noise in the Archeops data were used when com- 
puting the Archeops CMB angular power spectrum in 
(Bcnoit et al 2nn3all . 

5. Wavelet destriping 

Data from CMB experiments are commonly limited by 
intrinsic l//-type noise in the detectors. For experi- 
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Fig. 4. Top: Reconstructed time modulation function (j{t) 
of the Archeops noise in bolometer 217K04. Bottom: Time 
modulation function a{t) for a non-stationary simulation 
of the Archeops noise in bolometer 217K04. 

ments with a circular scanning strategy, like Archeops, 
WMAP and Planck, the l//-type noise shows up like 
stripes in the sky maps. The techniques used to remove 
those stripes are in general called destriping algo- 
rithms dEfstathiou, G 2005| 



Maino, D. et al. 20021 



Keihanen, E. et al 2004 



Poutanen, T. et al. 2004 



Revenu, B. et al. 20001 jSbarra, C. et al. 2003| | 



These 

are both based on the statistical properties of l//-type 
noise and on the fact that the noise, by contrast to the 
sky signal, is not coherently projected in the maps. 

5.1. Derivation of a wavelet destriping algorithm 

In the previous sections we have shown that the wavelet 
transform nearly diagonalizes the covariance matrix of 
l//-type noise. This allows us to design a destriping algo- 
rithm based on the wavelet transform. 

We assume that the time domain data for a typical 
CMB experiment can be written as 



dt = St 



(44) 



where St is the sky signal, and and n'^ are respectively 
a white and a correlated (color) noise component. We ex- 
press the projection of the time domain signal, St into a 
sky map, mp, via a pointing matrix Ppt such that 



rrLp = Ppt St 
Therefore we can rewrite equation lUljl as 



dt = Pt 



tpTflp 



(45) 



(46) 



The main purpose of a destriping algorithm is just to re- 
move or reduce significantly the contribution from the cor- 
related noise. This is obtained by maximizing the likeli- 
hood function, £, over the noise and the sky signal. 

We can characterize the correlated noise n'^ using its 
wavelet transform 



(47) 



for which the wavelet coefficients are Gaussian distributed 
and quasi-uncorrelated. We can then write the Hkelihood 
function as follows 



C = P{d\m)=P{W.^\m)P{Wnc\m) 



exp 



(48) 
(49) 



where and C„c are the covariance matrices in wavelet 
space of the correlated and white noise respectively 



(50) 
(51) 



Fig. 5. Baseline reconstruction on simulated Archeops 
data using a wavelet based destriping algorithm. In black 
we trace the simulated Archeops TOD. We overplot in 
blue, green and red the reconstructed basehne for the first, 
second and third step of the algorithm respectively. 



Notice that for our assumptions both matrices are diago- 
nal. 

By taking the log of the likelihood function (I4H|I we 
obtain 



A w^iu w I nc ^ nc ^ 



(52) 
(53) 

. . . (d - W^Wnc ~ Vm) + Wl,C-^Wnc (54) 

whose minimization with respect to m and Wnc leads to 
the following system of equations 

m = {VV)~^V'^{d-W'^Wnc) (55) 
= -C-^W{I -V{VV)-^V'^){d-W^Wnc) + 

■ ■■C-lWnc (57) 



which can be combined into 

(w(/ - rivvr^v^)w^ + Cy,c~,^)Wr. 
yv{i ~v{vv)-^v^)d 



(58) 



The solution to this system of equations can be found 
using iterative methods and in particular we use a conju- 
gate gradient method l |Press, W.H. et al. 1996| |. 

5.2. Practical implementation 

For 1 / /-type noise the variance of the wavelet coefficients 
at level j is given by 



(59) 



where cr^^ is the white noise variance and cr^ 2^°' is the 
variance of the correlated part with a being the 1// ex- 
ponent as described in equation 1271 

When the noise dominates the sky signal, we can ex- 
tract the noise parameters Gw, <J and a directly from 
data itself via Monte Carlo Markov Chain methods 
l lWada S. and Ito N., 20001 . 

When the signal contribution can not be neglected, the 
data wavelet variance is given by 



rr 2,2 



(60) 



Gb J and aw ,j can be iteratively evaluated from data 
using a template of the sky signal which is improved at 
each iteration. From the wavelet transforms of the tem- 
plate and the data we obtain 



~sky ,j 



(61) 
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Fig. 6. Top panel: Destriped map of the simulated 
Archeops data at 545 GHz using the wavelet destriping al- 
gorithm. Bottom panel: Residual stripes on the destriped 
map above. 

5.3. Application to simulated data 

The wavelet destriping was applied to simulated Archeops 
TOD at 545 GHz with 3 x 2^1 samples. We considered 
two main components in data, Galactic dust emission and 
l//-type noise. For the former we used a template of the 
Galactic dust emission scaled down to 545 GHz provided 
bv IFinkbeiner et al 19991 The l//-type noise properties 
were deduced from the Archeops data. 

We proceeded in three main steps. First of all, we esti- 
mated the lowest frequency components of the data from 
its wavelet decomposition and subtract it. Secondly, we 
improved this estimation by solving equation [^Hl itera- 
tively starting with largest scales and adding smaller and 
smaller scales progressively. Finally, we computed an ap- 
proximation to the signal, St, by thresholding the previous 
destriped map and deprojecting into the time domain. We 
then, step two was applied to the residuals, dt — st. 

Figure shows the simulated Archeops TOD at 
545 GHz. In blue, green and red we overplot the recon- 
structed noise low frequency components for step one, two 
and three respectively. We observe that we improve signif- 
icantly the estimated baseline in step three reducing the 
stripes in the maps as shown in figure El where we repre- 
sent the destriped map on the top panel and the residual 
stripes on the bottom panel. 

6. Time-frequency visualization and analysis of 
the Archeops TOD. 

In the previous sections we concentrated on the wavelet 
description of random Gaussian processes. We studied the 
decorrelation properties of the DWT and their application 
to the statistical analysis of CMB data sets and in par- 
ticular to the destriping of CMB maps. From the point 
of view of CMB data analysis this is of great interest but 
it is not all we can obtain from the wavelet transform. 
Actually, the most important property of wavelets is their 
simultaneous time and frequency localization. At this re- 
spect, wavelet analysis and in particular the DWPT is 
a fundamental tool for data visualization and character- 
ization. CMB time series dt are in general a linear com- 
bination of signal, both galactic and cosmological in ori- 
gin, of systematics effects like atmospheric contamination, 
and/or parasitic noises, and/or electromagnetic contami- 
nation, and/or glitches, and of random Gaussian noise 

aP^^sys^ + a'^'sysf" + . . . + nt ^ ' 

These components are different in nature and they present 
different time and frequency properties. For example, the 



Fig. 7. Top: Time-frequency analysis of the TOD of the 
217K04 Archeops bolometer using the DWPT. Bottom: 
Time-frequency analysis of the expected Galactic emission 
in the TOD of the 217K04 Archeops bolometer using the 
DWPT. See text for details. 

Fig. 8. Top: Time-frequency analysis of the TOD of the 
217T04 Archeops bolometer using the DWPT. Bottom: 
Time-frequency analysis of the expected Galactic emission 
in the TOD of the 217T04 Archeops bolometer using the 
DWPT. See text for details. 

cosmological and galactic signals only depend on the po- 
sition of the sky observed which is fully given by the in- 
strument pointing. By contrast, the atmospheric signal 
depends both on the sky position and on the time evo- 
lution of the atmospheric conditions. Further, other sys- 
tematic effects vary mainly with time and in general are 
related to physical phenomena in the detector surround- 
ings. Finally, the noise component is intrinsic to each de- 
tector and it is often of l//-type, Gaussian and locally- 
stationary. As the cosmological signal, the CMB emission, 
is in the time domain much weaker than all other com- 
ponents, these have to be identified, characterized and re- 
moved with high degree of precision before the extraction 
of the CMB anisotropies power spectrum. 

6.1. Time-frequency analysis of the Archeops data 

In the following we present the most relevant issues in the 
wavelet analysis of the Archeops data set. A more detailed 
description of the Archeops data and its processing can 
be found in ( [Macias-Perez, J.F. et al. 2005| |. To simplify 
further discussions we will assume only four main compo- 
nents in the Archeops timelines: CMB emission, galactic 
and atmospheric emissions, non identified systematics and 
Gaussian noise. 

The top panel of figure El shows the DWPT time- 
frequency representation of a typical Archeops detector 
signal. We observe mainly two components which dom- 
inate at high and low frequencies respectively. The low 
frequency signal which varies significantly with time is 
mainly given by galactic emission as shown by the bottom 
panel of the same figure where we represent the expected 
galactic emission for that detector. The circular scanning 
strategy of Archeops is such that during most of the ob- 
servation time we cross the Galactic plane perpendicularly 
and therefore the galactic emission shows up as a spike in 
time. By contrast in the last hours of fiight, the scanning 
direction is collinear to the Galactic plane and then, the 
galactic emission becomes broader in the time direction 
although we can still observe few frequency spikes related 
to intense and compact regions of the Galactic plane. The 
nearly time-fixed structures observed at low frequency al- 
though they look very much as 1// noise, are due to sys- 
tematics mainly dominated by atmospheric emission. 
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Fig. 9. Top: Average power spectrum of the TOD of the 
217K04 Archeops bolometer as a function of frequency 
computed from its DWPT Bottom: Time evolution of 
the power spectrum of the TOD of the 217K04 Archeops 
bolometer computed from its DWPT. 



The high frequency signal corresponds to the detector 
white noise whose power is artificially raised up with 
increasing frequency when deconvolving from the bolome- 
ter time constant. We observe that the noise smoothly 
decreases with time showing a clear non-stationary 
behavior. Indeed, as time went on during the fiight the 
bolometer temperature decreased therefore reducing the 
bolometer intrinsic noise. A detailed wavelet description 
and modehng of the non-stationarity of the Archeops 
noise is presented in the following subsection. 

The top panel of figure|Hlshows the DWPT of the worst 
Archeops bolometer (217T04) both in terms of noise and 
systematics. As above we plot in the bottom panel the 
DWPT of the expected Galactic signal for this bolometer. 
We observe that the DWPT of 217T04 presents as above 
at low frequencies the Galactic and atmospheric compo- 
nents and at high frequency a time varying noise. However 
it is dominated by a highly time varying signal in the fre- 
quency range 25 to 35 Hz. This signal is not observed 
in the data presented above and can be clearly identified 
as the residuals from unknown systematic effects. We can 
therefore, just by visual inspection, characterize the qual- 
ity of the Archeops bolometers data and identify which of 
them are badly affected by systematics. We have used this 
technique to identify the best, in terms of sensitivity and 
low level of systematics, Archeops bolometers which were 
used to construct the Archeops sky maps. 

6.2. Modeling of the non-stationarity of the Archeops 
noise 

As discussed in the previous section, the Archeops cryo- 
stat radically increased temperature when taking off and 
then cooled down slowly to achieve a nominal temper- 
ature of 95 mK for the last ten hours of fiight. This 
produced a fair decrease of the noise power with time 
which is the cause of the non stationarity of the Archeops 
bolometer noise. To account for this non stationarity we 
modeled Archeops noise as a time modulated-stationary 
wavelet process as it can be considered as locally (piece- 
wise) stationary noise with slowly time varying power but 
for the first two hours of fiight. Following subsections 14.21 
and 14.41 we have computed the a{t) function for each of 
the Archeops bolometers. The top panel figure El shows 
mean wavelet power as a function of frequency for the 
the Archeops bolometer 217K04. The bottom panel shows 
the reconstructed (j(t) function for the same bolometer. 
From these two quantities we can produce non-stationary 
simulations of the Archeops noise as shown in figure 21 



Fig. 10. Top: Wavelet destriped map of the Archeops 
545 GHz after one iteration of the algorithm described in 
section bottom; Wavelet destriped map of the Archeops 
545 GHz after two iteration of the algorithm. 

Fig. 11. Performance of the wavelet filtering in the 
353 GHz Archeops data. In orange we trace the recon- 
structed data baseline using a fit to a base of Gabor atoms. 
In blue we represent the reconstructed baseline using the 
wavelet detrending algorithm. 



These non stationary simulations of the Archeops bolome- 
ter noise were used to determine the angular power spec- 
trum of the noise in the Archeops maps which were needed 
for computing the Archeops CMB angular power spectrum 
in HBenoit etal 2nn3ajl via a MASTER-like algorithm. 

6.3. Wavelet destriping of the Archeops data 

We have applied the wavelet destriping algorithm pre- 
sented in section to the Archeops data at 545 GHz. 
When working with the simulated Archeops data we have 
only considered two main components. Galaxy emission 
and l//-type noise. As shown above, for the Archeops 
data we also need to consider the atmospheric emission 
at low frequencies which mimic a l//-type but it is 
not Gaussian distributed. To overcome this problem, we 
have first estimated and removed a very low frequency 
baseline in the data. After that, we have applied two 
iterations of the wavelet destriping algorithm as described 
in section [3 The results are presented in figure EH The 
top figure corresponds to the destriped map after one 
iteration of the algorithm. Residual stripes dominate 
the right top corner of the map. These are significantly 
reduced by the second iteration at shown in the bottom 
plot. It is important to notice that destriping allows us 
to recover the diffuse Galactic emission on the Gemini 
region (middle left of the figure) . 

The wavelet destriping algorithm presented here is 
mainly based on the minimization of the variance per pixel 
in the final map. As discussed in [Bourrachot, A. 2004| and 
[Macias-Perez, J.F. et al. 2005| this approach is not opti- 
mal for the Archeops data because of the atmospheric 
emission. A better approach is to minimize in the map 
the ratio between the variance perpendicular and paral- 

lei to the scanning direction, -|-. The generalization of 
the wavelet destriping algorithm to this latter approach 
is straight forward. As its implementation is much harder 
than when considering Gabor atoms as base functions, it 
has not being considered for the analysis of the Archeops 
data. However, despite of a simpler approach the qual- 
ity wavelet destriped maps can be compared to that of 
the final destriped maps used for the Archeops analysis 
dMacias-Perez, J.F. et al. 2005| |. 
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Fig. 12. Top: Fourier filtered map of the Ardieops 
353 GHz data Bottom: Wavelet filtered map of tlie of the 
Archeops 353 GHz data. See text for details. 

6.4. Wavelet detrending of the Archeops data 

As shown above, the Archeops destriped maps contain 
residual stripes which need to be removed before any 
scientific analysis of the data. In general these stripes 
are superposed to the sky signal of interest and therefore 
a careful filtering is needed. Further, as the Galactic 
signal in the Archeops maps is spike-like any filtering 
will produce ringing in the final maps. To overcome these 
problems, we have developed a wavelet based detrending 
algorithm. First of all, we mask the Galactic signal 
and obtain a first approximation of the low frequency 
components of the data by fitting them to a base of 
Gabor atoms. Then, we use the fitted data to interpolate 
over the Galactic mask. Finally, via a wavelet denoising 
algorithm we obtain the trends on the interpolated 
data which are then removed before map making. For 
denoising we use a wavelet thresholding algorithm limited 
to the few first smaller scales, typically up to j = 9 when 
working with G-miUions data samples. 

Figure [TTl shows the performance of the detrending al- 
gorithm in time domain. We plot in orange the interpo- 
lation to the data obtained by fitting Gabor atoms and 
in blue the baseline obtained via the wavelet algorithm. 
We clearly observe that the the wavelet baseline manages 
to follow much better the accidents on the data reducing 
both the residual stripes and ringing in the final maps. The 
bottom panel of figure El represent the wavelet filtered 
map at 353 GHz which can be compared to the Fourier 
filtered map in the top panel of the figure. In the wavelet 
filtered maps the residual stripes and ringing are much 
less. Further the diffuse structure of the Galactic emission 
is better preserved. Notice that the wavelet filtered maps 
at 353 GHz were used for the determination of the polar- 
ization of the diffuse Galactic dust emission with Archeops 
ijBenoit et al 20M\ . 

7. Conclusions 

The discrete wavelet transform (DWT) , maximum overlap 
discrete wavelet transform (MODWT) and the discrete 
wavelet packet transform (DWPT), because of their 
simultaneous time and frequency localization properties, 
are very important tools for the analysis of TOD from 
CMB experiments. They allow us to trace the evolution 
in time of the data power spectrum and to easily identify 
systematic effects on the data. 

The DWT permits a compact representation of 
Gaussian stationary noise. Indeed, it decorrelates 
Gaussian l//-type noise which is commonly present in 
the detector of CMB experiments. Within the wavelet 



description the covariance matrix of Gaussian l//-type 
noise is diagonal as it is the case in the Fourier space. 
This allows us to efficiently and accurately simulate 
Gaussian l//-type noise using wavelets. Moreover, the 
DWT transform permits a straight forward description of 
locally stationary Gaussian noise via the time evolution of 
the variance of the wavelet coefficient but with a diagonal 
covariance matrix. A particular case of this is the time 
modulated Gaussian noise for which the time evolution is 
common to all the wavelet scales. 

The above properties allows us to generalize the 
Fourier space algorithms for fast optimal map making 
and maximum likelihood determination of the CMB 
angular power spectrum to the wavelet space including 
both stationary and locally stationary Gaussian noise. At 
this respect, we have developed a wavelet based destriped 
algorithm which both in simulated and true Archeops 
data reduce significantly the level of stripes in the maps. 
Despite the fact that this algorithm is based on a simple 
approach of minimization of the variance per pixel, the 
results obtained can be compared to those from more 
precise destriping algorithms. 

Finally we have performed a full wavelet analysis of the 
Archeops data. The visuaHzation and careful study of the 
time-frequency space via the DWPT allows us to clearly 
identify systematics and to characterize the quality of the 
data. Further, we have proceeded to the careful wavelet 
modeling of the locally stationary noise in the Archeops 
data. This modehng allowed us to obtain via simulations 
the angular power spectrum of the Archeops noise needed 
for the estimation of the CMB angular power spectrum 
with Archeops IjBenoit et al 2nfl3a)l . We have applied the 
wavelet destriping algorithm to the Archeops data obtain- 
ing encouraging results. Finally, we have developed a de- 
trending algorithm based on a wavelet denoising of the 
data. This algorithm applied to the Archeops destriped 
data reduce significantly the residual stripes on the final 
maps and introduces very little ringing. The wavelet de- 
trended Archeops maps at 353 GHz were used for the de- 
termination of the polarized diffuse emission from Galactic 
dust l|Benoit et al 2(\M\ . 
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